Mon. Not. R. Astron. Soc. 000,[T](24](2009) Printed 15 November 2010 (MN KTbX style file v2.2) 



Application of XFaster power spectrum and likelihood estimator 
to Planck 

G. Rocha^'^, C. R. Contaldi^ J. R. Bond^ & K. M. Gorski^'^ 

^ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena CA 91109, U. S. A. 
^ Department of Physics, California Institute of Technology, Pasadena, 91125, U. S. A. 

^ Department of Physics, Imperial College, University of London, South Kensington Campus, London SW7 2AZ, U. K. 

* The Canadian Institute for Theoretical Astrophysics, CITA, University of Toronto, 60 St. George Street, Toronto, Ontario, M5S 3H8, Canada 
^ Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland 

27 October 2010 



ABSTRACT 

We develop the XFaster Cosmic Microwave Background (CMB) temperature and polar- 
ization anisotropy power spectrum and likelihood technique for the Planck CMB satellite 
mission. We give an overview of this estimator and its current implementation and present 
the results of applying this algorithm to simulated Planck data. We show that it can accu- 
rately extract the power spectrum of Planck data for the high-^ multipoles range. We compare 
the XFaster approximation for the likelihood to other high-i' likelihood approximations 
such as Gaussian and Offset Lognormal and a low-£ pixel-based likelihood. We show that the 
XFaster likelihood is not only accurate at high-£, but also performs well at moderately low 
multipoles. We also present results for cosmological parameter Markov Chain Monte Carlo 
estimation with the XFaster likelihood. As long as the low-^ polarization and temperature 
power are properly accounted for, e.g., by adding an adequate low-i likelihood ingredient, the 
input parameters are recovered to a high level of accuracy. 

Key words: Cosmology: observations - methods; data analysis - cosmic microwave back- 
ground 



1 INTRODUCTION 

Power spectrum estimation plays a crucial role in CMB data analysis. Primordial curvature fluctuations form a homogeneous, isotropic, and 
nearly Gaussian random field in most early universe scenarios, inflationary or otherwise. To the extent that fluctuations are Gaussian, the 
power spectrum describes their statistical properties fully. An immediate consequence is that the CMB temperature and polarization primary 
anisotropies, linearly responding to the primordial fluctuations, form an isotropic nearly Gaussian random field, characterized by their own 
angular power spectra. 

With the advent of large, high-quality data sets, especially that from the Planck mission ( [Planck Blue Book|j2005) ; |Tauber et aTlpOOgt ), 
we can measure the angular power spectrum accurately over a wide range of angular scales. Power spectrum estimation can be viewed as a 
significant data compression. For instance, 1 year of observations from a Planck High Frequency Instrument (HFI) channel produces roughly 
Ntod ~ 2 X 10^^ data samples. This reduces to Nmap ~ 5 x 10^ pixels via mapmaking and finally to Npse ~ 3 x 10'^ values in the power 
spectra. 

Planck is a full-sky experiment with beams ranging in size from 30' to 5', and with high resolution maps encompassing tens of millions 
of pixels. Direct extraction of science from the pixelized maps is computationally expensive and in fact unfeasible. Accurate estimation of 
the angular power spectrum for Planck enables the extraction of science with minimal loss of information. 

A number of approaches have been developed to estimate the angular power spectrum from CMB data (for a review see Efstat hioul 
( |2004^ ; [Ashdown et al.| ( (20To] )). Such estimators can be divided into three classes: codes accurate at large-angular scales (low multipoles i), 
including codes that evaluate or sample from the likelihood function directly; codes accurate at small angular scales (high-^) that characterize 
the statistics of an unbiased frequentist estimator for the power spectrum; and hybrid codes that can be aplied to both low and high-^. The 
first class comprises Maximum Likelihood Estimators (MLE) either in Fourier space (G6rski|( (l994a[|1997) ; |G6rski et al.H1994b[fl996[ l) or 
in Real space ( [Tegmark, Bunn| ( |1995^ ; |Hancock et al.| ( [1997[ >) such as MADspec (^orrill ( 1999); |Borrill et al.|p009) ). Quadratic Maximum 
Likelihood (QML) estimators jWright et alTl ( [T994) ; [Hamilton| ( [TW7l >; |Tegmark|p997^ ; |Bond" Ja^^^ & Knox| ( (2000l >; [Knox] (1999) ) such as 
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Bolpol fGruppuso e"tall ( |2009| l), Gibbs samplers such as Commander jjewell et alT] l |2004) ; [Enksen et al.| ( |2004^ ; |Wandelt et al.|f2004^ ), and 
Importance Samplers combined with a Copula-based approximation to the likelihood such as Teasing ([B enabed et al.| j2009|). The second 
class comprises quadratic Pseudo — Ce (PCL) or Master method codes ( [Hivon, Gorski et ar] ('2002 1; Wandelt, Hivon & Gorski ( 200TJ) such 
as Romaster, cRomaster jPolenta et al.|p004l l), Xpol jTristram et al.|p005a|b| l), and Crosspect ( Ashdow n et al., (,2010)), angular correlation 
function codes such as Spice ( |Szapudi et al.| ( |2001^ ) and Polspice ( [Chon et"aL] ( |2003^ ), and Quadratic Maximum likelihood (QML) codes 
such as XFaster ( [Contaldi et al.| ( |2010^ ; |Netterfield et al|p002l l; |Montroy et al| ( (2006l l and this paper, see Section [3TT] for more details). 
The third class consists of hybrid power spectrum estimators such as a QML estimator at low-^ combined with a PCL estimator at high-^ 
( |Efstathioul ( |2004| l) and, to some extent, the XFaster method alone. 

XFaster ( [Contaldi et al.|p010| l) was first developed to give rapid and accurate power spectra determinations from bolometer data for 
the Boomerang long-duration balloon experiments, first for total anisotropy ( [Netterfield et al[ ( (2002| l) and then for polarization ( [Montroy et al[ 
( [2006^ ). XFaster is a quadratic maximum likelihood estimator formulated in the isotropic, diagonal approximation of the Master method 
( [Hivon, Gorski et al.[(2002l >). The noise becomes a diagonalised Monte-Carlo-estimated bias and the signal is summed into bands to reduce 
correlations induced by sky cuts. In this sense XFaster is an extension of the traditional Master estimators where the pseudo-C^ quantity 
is replaced by the quadratic MLE expression and uncertainties are given by the Fisher matrix. 

This method has been compared with other high-i" codes such as Polspice, Romaster, Xpol, CrossSpec within Planck working group 
C'l Temperature and Polarization, CTP. A detailed account of this comparison will be given in a paper by the Planck CTP working group 
( [Ashdown et al.[ ( |20I0^ ). A full and detailed account of XFaster as a standalone method will be given elsewhere ( [Contaldi et al.[pblO^ ). 
Here we give an overview of the method, but our main goal is to show its adequacy to extract the power spectrum from Planck data. 

An interesting feature of the method is that it provides a natural expression for the likelihood based on the assumption that the cut-sky 
harmonic coefficients aim follow the same distribution as those of the full-sky harmonics. We compare our approximate likelihood to the 
exact full-sky likelihood (the inverse Wishart distribution) and to the pixel based likelihood (i.e., multivariate Gaussian of the pixel's /, Q, 
and U Stokes parameters). We show that XFaster agrees well with the exact likelihoods at moderate low-£ multipoles as well. 

Using the XFaster power spectrum and likelihood estimator, we show how to go straight from the map to parameters, bypassing the 
band power spectrum estimation step. Alternatively, we show how to use the band power spectrum estimated with XFaster in combination 
with any likelihood approximation to estimate parameters. In particular, in [Rocha et al.[ ( [2010[ l we compare parameters estimated with 
XFaster and the Offset Lognormal Bandpower likelihood to those obtained with the XFaster likelihood. From our analysis we conclude 
that as long as the low-£ polarization is properly accounted for (by adding an adequate low-i" likelihood ingredient), we recover the input 
parameters accurately. 

This paper is organized as follows: Section|2]describes the map and the Monte Carlo simulations for two phases of increasing complexity 
studied in the CTP working group; Section [3] gives an overview of the XFaster power spectrum and likelihood estimator, including 
the estimation of kernels, transfer or filter functions, and window functions; Section |4] shows the results of applying XFaster to Planck 
simulations in several different ways. It includes the impact of beam asymmetries on power spectrum and cosmological parameter estimation, 
comparison of the XFaster likelihood to other likelihood aproximations at high-^ and to pixel-based likelihood at low-l, and cosmological 
parameter estimation. Section|5]gives conclusions. 



2 SIMULATIONS 

The Planck satellite ( [Planck Blue Book|j2005| >; |Tauber et alH2009| ) is a full-sky experiment with beams ranging in size from 30' to 5'. The 
Low Frequency Instrument (LFI) covers 30, 44, and 70 GHz; the High Frequency Instrument (HFI) covers 100, 143, 216, 353, 545, and 
857 GHz. From the second Lagrangian point of the Earth-Sun system (L2) Planck scans nearly great circles on the sky, covering the full sky 
twice over the course of a year ( [Pupae & Tauber ] ( [2005| >). Planck spins at I rpm around an axis that is repointed roughly 30 times per day 
along a cycloidal path, with the spin axis moving in a 7? 5 circle around the anti-Sun direction with a period of six months. This ensures that 
all feeds cover the ecliptic pole regions fully. We also include small perturbations to the pointing, with spin axis nutation and variations in 
the satellite spin rate. For the analysis presented here we consider the 70 GHz LFI channel. 

The simulations used in this work include CMB and realistic detector noise only, and are specified by the scanning strategy (as described 
above), telescope beams, and detector properties. To mimic a more sensitive combination of channels, the white noise level was taken to be 
lower than that expected for the real 70 GHz channel. We used a single observed map containing CMB and noise as well as Monte Carlo 
simulations of signal and noise. Technical details of the simulations are given in [Ashdown et"^ ( [20T0l l. 

We considered all twelve detectors of the 70 GHz LFI channel. The beams of the detectors have FWHMs of I3'~14', so the maps were 
made with Nsidc ~ 1024, corresponding to a pixel size of 3.'4. Two sets of maps were provided, one I2-detector map to be used in the 
auto-spectrum mode, and three 4-detector maps to be used in the cross-spectrum mode. 

The input sky signal used to generate the observed map was the CMB map derived from the WMAP 1-year data used in a previous CTP 
map-making exercise jAshdown et al. ( 2009 2010jl). It is derived from the Planck CMB reference sky available in[^ Hence the large scale 



structure of the observed map is a WMAP constrained realization. The angular power spectrum of the aim is plotted in Figure[T| 



^ http://www.sissa.it/ planck/reference-sky/CMB/alms/alm-cmb-reference-template-microKthermodynamic-nside2048.fits 
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Simulations in four steps of increasing complexity were used. For historical reasons we refer to these steps as Phases la, lb, 2a, and 2b. 

Phase la — Data were simulated with symmetric beams and isotropic white noise. The power spectrum was estimated from the full sky. 

Phase lb — Data were simulated with symmetric beams and anisotropic white noise determined by the scan strategy. A sky cut was applied 
in the calculation of the power spectrum to mimic the effects of removing the galactic plane from observed data. 

For Phases la and lb the input CMB sky was convolved with a Gaussian beam of 14' FWHM, and the noise realizations were generated 
assuming an RMS of 69.28 /iK per pixel in temperature and am RMS of 97.97 /iK per pixel in the Q and U polarization components. 

The observed maps were generated in the pixel domain. Monte Carlo signal simulations were generated from the WMAP 1-year best-fit 
ACDM power spectrunj^ plotted in Figure [T| (blue solid line). One hundred Monte Carlo realizations were generated directly in the pixel 
domain for both the signal and noise using HEALPix tools ( Gorski et al.| ( |2005,) ) and our own simulator 

Phase 2a — Data were simulated with both correlated 1// and anisotropic white noise. Symmetric beams were assumed, all Gaussian with 
FWHM 14'. Noise effects induced by temperature fluctuations of the 20-K hydrogen sorption cooler were also included. 

Phase 2b — Data were simulated with both correlated 1// and anisotropic white noise. Asymmetric beams were used, specifically, elliptical 
Gaussians fit to the central parts of realistic beams calculated by a full diffraction code for the Planck optical system. For the twelve beams, 
the geometric mean of the major and minor axis FWHMs ranged from 12.'43 to 13.'03. Major axis to minor axis ratios varies from 1.22 to 
1.26. 

For Phases 2a and 2b the white noise per sample was 2025.8 /iK; the 1// noise power spectrum had a knee frequency of 0.05 Hz and a slope 
of -1.7. 

The observed maps were made from time-ordered data (TOD) using the destriping mapmaker Springtide ( |Poutanen| ^2005^ ; | Ash- 1 
[down et al.| ( [2007a|b||2009) ; |Ashdown| |2009c)). The TOD were generated using modules of the Planck simulator pipeline. Levels ( |Rei-| 
|necke et al.| ( |2005^ ). Where a sky cut was applied in the analysis of the maps, the cut was made at the boundary where the total intensity 
of the diffuse foregrounds and point sources exceeded twice the CMB sigma. Masks for missing pixels due to the scanning strategy, if any, 
were also considered. Figure |2] shows the observed map for Phase2b using all twelve detectors, with the mask for galaxy plus missing pixels 
applied. 

As with Phases la and lb, 100 Monte Carlo signal simulations were generated from the first year WMAP-l-CBI-l-ACBAR best fit ACDM 
power spectrum^ with BB mode power set to zero (see Figure [!}. For the symmetric beam case only the noise TODS were generated, while 
the signal was simulated in the map domain. For the asymmetric beam case both signal and noise simulations were generated in the time 
domain. 

As mentioned above the large scale structure of the observed map is derived from real observations, i.e., a WMAP constrained re- 
alization, hence it is not necessarily consistent with the best-fit spectrum at low multipoles. This discrepancy will become evident later 
when comparing the power spectrum estimated from the observed map with the best-fit theoretical spectrum, as well as when compar- 
ing the cosmological parameters estimated with XFaster power spectrum and likelihood and the theoretical best fit parameters. As the 
Monte Carlo simulations are realizations of the WMAP 1-year best-fit ACDM power spectrum for Phases la and lb, and the first year 
WMAP-l-CBI-l-ACBAR best fit ACDM power spectrum for Phases 2a and 2b, such discrepancy is no longer present. Parameters estimated 
from these Monte Carlo simulations maps are now close to the WMAP best fit parameters. 

The choice of 70 GHz for the simulations was driven by practical matters of computational resources having to do with the size of 
the TOD, the number of pixels in the maps, and the number of multipoles that had to be calculated. The HFI channels have higher angular 
resolution and sensitivity, and will extend to smaller angular scales with reduced error bars. Increases in computational capability over time 
make it possible now to generate thousands of Monte Carlo simulations at the higher frequencies as well. Results will be presented in a future 
publication |Rocha et aLl ( |2010| l; [Ashdown et aLl ( |2010^ . 



3 XFAS TER POWER SPECTRUM AND LIKELIHOOD ESTIMATOR 
3.1 XF a s t e r Power spectrum estimator 

XFaster ( jContaldi et al.|p010^ ; |Netterfield et~al|j2002) ; |Montroy et al|j2006^ ) is an iterative, maximum likelihood, quadratic band power 
estimator ( |Bond Jaffe & Knox|(2000l >; |Knox| ( [T999l l; |Hamiltonl ( [TW7) ; [Tegma^^^ & de Oliveira-Costa| ( |200T] l) based on a 

diagonal approximation to the quadratic Fisher matrix estimator. It is a QML estimator formulated in the isotropic, diagonal approximation 
of the Master method ( |Hivon, Gorski et al.|f2002t ). 

It is common to expand the pixel temperature fluctuations (Stokes /), r(n), on the celestial sphere in terms of spherical harmonic 
functions, Y^m, as 

T(n) = ^a,„y<™(n), (1) 



available in http://lambda.gsfc.nasa.gov/product/map/drlAcdm.cfin 
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Figure 1. Angular power spectrum, Ci, of the Planck CMB reference sky, obtained from the reference sky spherical harmonics, a^^, used to generate the 
signal of the 70 GHz observed map (green solid line); used to generate the Monte Carlo simulations for Phase 1 (blue solid line); Ci used to generate the 
Monte Carlo simulations for Phase 2 (black solid line). 




Figure 2. 70 GHz (/, Q, U) Stokes parameters of the observed map (from left to right hand side) for Phase 2b, generated with all detectors, convolved with 
the asymmetric beam and with mask (for galaxy plus missing pixels) applied. 



with coefficients aim ■ 

The maximum likelihood estimator (MLE) is based on a Gaussian assumption for the hkelihood of the observed data (e.g., the pixel 
temperature, T, or its spherical harmonic transform a^m): 

1 



L(d|p) 



(27r)'v/21Cli/2 



exp 



c 



(2) 



where C is the covariance of the data, and p is the set of model parameters. C(p) = S(p) + N, where S is the sky signal and N is the noise. 
For single dish, full sky observations an isotropic signal is diagonal in the spherical harmonic space and can be described by an m-averaged 
power spectrum Ci on each multipole, i.e., Sirn.i'm' = 5ui Smm'Ci. The noise is generally not diagonal. 

All high-^ codes assume the data to be Gaussian distributed, however, they differ from XFaster on the devised unbiased frequentist 
power spectrum estimator. All algorithms form quadratic functions of the data. Pseudo-C^ (PCL) codes estimate Ci = X/m=-f ^"4" 
using fast spherical transforms. Example of such algorithms are all Master type codes (Hivon, Gorski et aL]p002| l;[Wandelt, Hivon & Gorski | 
((20011] such as Romaster, cRomaster ( Polenta et al.|^2004^ ), Xpol ( [Tristram et"ar] ( |2005a|bt ), and Crosspect ( Ashdown et aL| ( |2010| l). Others 
calculate the angular correlation function C{9) = ^ X/i>2(^^ ~'~ l)C'<A(cos0) (where Pe{cos9) is the Legendre polynomial and S is an 
angular separation on the sky) using fast evaluation of the 2-point correlation function such as [Spice ( [Szapudi et al.|pOOT) ), Polspice jChon] 
|et al. (20031], though this is achieved using fast spherical transforms. XFaster ( |Contaldi et al.|(2010[ >; |NetterfieId et al| ( |2002[ l; |Montroyl 
|etalj(^2006y), instead uses the Quadratic Maximum likelihood (QML) expression as derived below. 



It can be shown that the(^, m) space maximum likelihood solution for the power spectrum is given by (Bond Jaffe & Knox (2000 
|Hamilton| ( [T997l >; |TegmaA|p997) ; [Tbgmark & de Oliveira-Costa| ( [200T) ): 



XF a s t e r applied to Planck 5 



(3) 



where = aem'^i^m' '^^e quadratic in the coefficients of the expansion of the observed map and -^[^is the Fisher information 

matrix for the C'e parameters (= curvature matrix in the ensemble average limit), given by: 



= \ Tr 



as 1 as 1 

dCi dCe, 



(4) 



An iterative scheme can be employed to reach the maximum likelihood estimate for the Cf. start with an initial guess; compute evalu- 
ate Eq[3] However, matrix operations become prohibitive for dimensions larger than a few thousand. 

To circumvent this problem XFaster recasts the estimator in the isotropic, diagonal approximations of the Master methods jHivon,| 
|G6rski et al.H20d2) ; |Wandelt, Hivon & G6rskiH200l"| l) simplifying the calculations above. In this case the noise becomes a diagonalised 
Monte Carlo estimated bias and the signal is summed into bands to average down the correlations induced by any reduced sky coverage. In 
this case for a single mode, say temperature alone, the covariance of the observed cut-sky modes is approximated by 

Clm,e'm' ~ Su'Smm'{Ce + (Ne)), (5) 

where Ce is the cut-sky model power spectrum. In our case the cut-sky power spectrum is parameterized through a set of deviations qi from 
a template full-sky 'shape' spectrum , 

I' 

where Kui is the coupling matrix due to the cut sky observations (see section [3.1.1[ l, Fi is a transfer or filter function accounting for the 
effect of pre-filtering the data both in time and spatial domain (see section [3.1.2^ , and Bi expresses the effect of a finite beam. For the case 
where the spectrum is parametrized in bands we consider band power deviations qt 

Ci = Y^ qbCZ = H I] Ku'BlF,,C§Xbie'), (7) 
6 be' 

where XbW is a binning function. Assume for simplicity flat binning with Xb{f-) = 1 within the band and zero outside. The ML solution for 
the qb is 

.b - ^ E^- E(2^ + '^'wfkr^^''^ ^ 

where isotropy reduces the trace as Tr — >■ 9 describes the effective degrees of freedom in the maps (which may be reduced 

by additional weighting of the modes such as filtering or pixel weighting), and is related to the moments of the pixel weighting and the sky 
coverage, see |Hivon, Gorski et al.|p002^ , and it can be further impacted by the binning of the power spectrum: 

g = fskyAi— where fskyW^ = / W*(n)dn (9) 

where VK(n) is the window or mask applied to the data, fskyWi is the i-th moment of the arbitrary weighting scheme, and is the width of 
the multipole bins. 

The expression for the Fisher matrix is now given by 

TbV = \y{2l+l)g .^''^'y , (10) 

For polarization sensitive observations the data include the /, Q, and U Stokes parameters. As mentioned before, the / map is expanded 
in terms of spherical harmonics while the Q and U maps are expanded in spin-2 spherical harmonics, 2Yim, to obtain E and B (grad-type or 
curl-type) polarization coefficients: 

± iaf„)±2Y<m(n) (11) 
There are six spectra representing the six independent elements of the 3x3 covariance matrix of the (aj^, af^,af^) vector: 

cr = ^grc'f^^+ivr (12) 



^ It is traditional in the literature to use to represent the Fisher matrix and F its ensemble average, however to avoid confusion with Fi that denotes the 
filter or transfer function we resort to F for the ensemble average of the Fisher matrix. 



6 G. Rocha, C. R. Contaldi, J. R. Bond & K. M. Gorski 



c 



c 



EE p;(S)EE gg ~{S)BB ZrEE 

^ +qt -^C^i, ) + Ni 



El EE 7^ 

b 
b 

J-^TE^ipTE^^TE 
b 

b 



B ^,iS)BB gg ~(S)EE ~gg 

bi +1b -Cm ) + Ni 



(13) 
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The template shape matrices are defined using the various coupling kernels for the different polarization types. The transfer functions 
are distinct for each polarization type: 



.p^(S)EE 
,p^(S)BB 
p^(S)TE 



(S)TB 



c, 



{S)EB 
bi 



= J^Ku'FrBlclf^^Xbil), 
I' 

= J2^K„,F,'^^Bicf^^Xb{e), 
I' 

= J2^K,„FrB!,C^P^^Xb{l), 
I' 

= Y,.K,„FrBlcf^^xb{e), 

I' 

= J2>^K,,,Fi;'^BlC^P^^Xb{l), 
I' 

= - -K,,,)F^^Blc\f^^''xb{t). 



(18) 
(19) 
(20) 
(21) 
(22) 
(23) 



For simplicity the beam Bi is assumed independent of polarization. (In principle, it could also be treated distinctly for each of the T, 
E, and B modes.) The mask coupling kernels, Kui , and the two additional polarization mask coupling kernels, ±Ku' , xKui are defined 
in section [J. 1.11 

Extending the above formalism to polarization, the XFaster estimator takes a matricial form, implemented trivially since the matrix 
C is now block diagonal: 

C — >■ diagCDe^-^ , Df^.^+i, . . . , D^^^^), where each multipole's covariance is a 3 x 3 matrix: 



D, = I cr cf^ 



c 



C 



C, 
C, 



I 

EB 



(24) 



Similarly its inverse is a block diagonal of the inverses of matrices and therefore simple to compute. The noise covariance matrix is also 
of this form, the iV/^^ in each block diagonal is obtained by noise only Monte Carlo simulations. The band power deviations, qb, take now 
the following form: 



lb 



(25) 



and the Fisher matrix is now given by: 



^bb' = ^^(2^+l)»Tr 



D, 



(26) 



where the band index, b, spans bands in all polarization types. The derivatives of the signal matrices with respect to the deviations qt are 
given by: 



dS 

dqb 



C, 



(S)TT 

be 







as 

dqb 



C, 



cip-- 

(S)TE 




(27) 
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I (28) 

^(S)BB 



^U::TB = | \ , ^\b::EB = \ C',^'^^ \. (29) 



MS)EB 





These derivatives include contributions from both (+) and (— ) kernels in the case EE and BB due to the geometrical leakage. 

This estimator makes use of the Monte Carlo pseudo-Cf formalism of Master methods to estimate noise bias and linear filter functions. 
It requires noise only Monte Carlo simulations to estimate the noise bias and signal only Monte Carlo simulations to estimate the filter 
function. Contrary to the conventional pseudo-C^ based methods it does not requires signal+noise Monte Carlo simulations to estimate the 
uncertainty (variance) of the band powers, which is given by the Fisher matrix instead, a byproduct of the method. The Fisher matrix is 
computed self-consistently and runs over all band powers and polarizations. It takes into account all the correlations in the approximation 
used (coupling kernel and diagonal noise bias) and can therefore be used to compute, self-consistently, all ancillary information required in 
the estimation process, correlations, window functions, etc. In Figure[3]we plot the inverse of the fisher matrix for phase2, symmetric beam 
case. Note that the Fisher matrix is not diagonal. 

Furthermore XFaster can estimate both auto-spectra and cross-spectra jointly, using the full covariance of the aim's, via a multiple- 
map analysis. 

As mentioned above, the noise is generally not diagonal. For Planck, the noise is white to good approximation at small angular scales; 
however, at large angular scales instrumental characteristics such as 1// noise and thermal fluctuations combine with the scan strategy to 
produce significant off-diagonal correlations in the noise. Therefore, the XFaster approximation is not optimal at low-£ multipole range. 
We show in Sections [4.2| and [43] that XFaster is a very good approximation for £ > 30. We have no intention of using XFaster for Planck 
at low-£. Instead, we will combine one of the codes adequate at \ow-l (listed in Section[TJ with XFaster (or another high-^ estimator) into 
a hybrid estimator of the power spectrum that covers the entire multipole range. 

Previous CMB experiments have customarily binned power spectra in multipole bands. The main reason for this is that it "enhances" 
the signal-to-noise ratio of the spectra. It also averages down correlations due to the reduced sky coverage at the same time as the coupling 
matrices due to the cut-sky attempt to correct the cut-sky effect. The XFaster power spectrum can be computed multipole by multipole, 
i.e., for each £, or in multipole bands. The band power spectra are given in Section [4T| To estimate cosmological parameters we can use band 
power spectra with any high-i' likelihood approximation, e.g., with the Offfset Lognormal Bandpower likelihood presented in |Rocha et al.| 
( |2010^ . However, as the XFaster likelihood is estimated multipole by multipole we can bypass the band power spectrum estimation step 
and estimate parameters directly from the maps (via its raw pseudo-Cf). Slices of the XFaster likelihood and parameter constraints are 
given in Section|42]and Section [43] 

As kernels and transfer (filter) functions are an important ingredient in power spectrum estimation we describe next how they are 
computed within XFaster approach. 



3.1.1 Kernels 

The effect of masking the sky is to produce a power spectrum that is a linear combination of the full sky power spectrum multipoles on the 
sky. The coupling matrix due to the cut sky observations, Kui , encodes this effect, it only depends on the geometry of the mask or window 
and is easily computable. 

Considering a window function W{n), and ignoring the effects of beam convolution and filtering effects due to any pre-processing of the 
timelines, the ensemble averages of the cut-sky (C'e) and the full-sky angular power spectrum (Ct) can be related by: 

(Cf) =^ii-«.{C,0, (30) 
e' 

with coupling matrix, Kui given by: 

= ?^^^J(A/, f; 0,0,0) (31) 

r 

f £ £' £" \ 

where J(/, Z', ; 0, 0, 0) = I I is the 3j symbol, and is the power spectrum of the window function VK(n), that is 







W; = {2£ + l)Wt with: 



>V^ = ^^^|m„|' and Wim = j dnWih)Y,*m{A) (32) 
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Figure 3. Logarithm of the absolute value of the 'normalized' (ie set to 1 at the maximum value) inverse of the Fisher matrix (covariance matrix) for phase2, 
symmetric beam case. TT, EE, BB and TE modes are displayed sequentially from bottom left-hand side corner to the upper right-hand side corner along the 
diagonal. 



Wo = Wo' = 47r/iyW? and = ^^{21 + ^ A-k f,kvW2 (33) 
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Extending the above to polarized data we consider the additional polarization mask coupling kernels defined as follows: 



xA-«. = ;(2L + i)w^[ ; ^ : M ; ; : i(i + (-ir'--) (36) 

These kernels account for the leakage of power between E and B modes induced by the usage of the full-sky ±2y«m(n) basis on a 
cut-sky. 

In Figure|4]we plot the masks for temperature and polarization used in Phase 2 (see section|2|. 
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cb6_i_gal_iicst.fits: UNKNOWNl 



cb6_p_gal^cst,fits: UKKKOWNl 




Figure 4. Mask for Temperature (left hand side) and polarization (right hand side) used in Phase 2. The /gky 
Polarization. 



0.85 for Temperature and f^y^y ~ 0.73 for 



3. 1.2 Transfer or Filter functions 

To compute Fi we start with Eq.[7]for the signal only Monte Carlo simulations, and as Fi is assumed to be smooth over each multipole bin 
we move Fi out of the summation for each bin b, to get: 



Ct = Yl 1"^^' = '^'^^ Ku'Blc!,Xb{i') (37) 

b be' 

and proceed with the iterative scheme as one would to estimate qt but now we estimate the transfer function, Fb instead. This is achieved by 
fixing qt = 1, varying Fb and considering the signal only Monte Carlo simulations. 



Fb 



(Si) 



(38) 



where (Se) is the average of the signal only Monte Carlo simulations. 
Extending to polarization we have: 



(39) 



3.2 XFaster likelihood estimator 

A very attractive feature of the XFaster power spectrum estimator is that is naturally provides a likelihood, i.e., the Probability of the 
observed cut sky data given the model. In the XFaster approximation, considering only one mode (say Temperature alone), the likelihood 
takes the following form, up to a constant (where A means A estimated on the cut-sky): 



Ce 



{Ce + (Ne)) 



+ \n{Ce + {Ne)) 



where Ce is the cut-sky model power spectrum given by Eq.|7]in Section 
consider band power deviations qb- 

Extending to Temperature and Polarization we have: 



3.1 



InL : 
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(40) 



for the case where the spectrum is parametrized in bands we 



(41) 



3.1 



Where Nf and are given in section 

An interesting point to note is that XF a s t e r likelihood follows intuitively from the usual full-sky ideal case exact likelihood (an inverse 
Gamma distribution for temperature alone and an inverse Wishart distribution for temperature + polarization (see e.g., |Rocha et al.| ( [2010| l). 

Here we use one dimensional slices as an approximation to investigate the non-Gaussianity of the likelihood. One samples in each 
deviation qb direction individually around the maximum likelihood solution q^. This approximation is adequate if the band powers are 
not heavily correlated. Note that the likelihood slices are estimated along the bands and not along each I, and hence will be affected by the 
binning procedure. To compare XFaster likelihood to other approximations we make use of slices computed along the bandpower spectrum 
deviations, qb- Such likelihood slices for the 70 GHz observed map are plotted in Section [4~2| 
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When estimating parameters witii XFaster likelihood, by default estimated multipole by multipole, we do not make use of the band- 
power spectra. It is in this sense that XFaster likelihood allows to go straight from the maps to parameters bypassing the band power 
spectrum step. It only requires the raw pseudo-C^ of the observations plus the kernel and transfer function to relate the cut-sky pseudo-C^ to 
the full-sky Ct- 



3.2.1 Window functions 

To compare the theoretical power spectrum to the observed power spectrum and to estimate parameters, we must construct an operator for 
obtaining theoretical bandpowers from model power spectra Cj . Following 

AM = E 



Bond Jaffe & Knox 



2000 1 we define a logarithmic integral 



+ ^ f 



(42) 



which is used to calculate the expectation values for the deviations g(,(when a shape model, Cf is considered), or bandpowers Cb (when Cf 
is assumed to be flat). 

1 [W^Ct] 



where is the band power window function, andC^^^ = £{£ + l)Cf/2TT. 



We define normalized window functions to be 
= 1. 



By taking the ensemble average limit of Eq. ([8b and using the fact that 



{(Cf^ - TV,)) 
we obtain 
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Extending to polarization: 



Wt = 



4-7r 



(21 +■ 



I' 



where Wt 



D 



-1 as_ 



D7\ and gives the cut-sky response to the individual full-sky multipoles: 



KinFf^Sf xKiifFf^ Bf xKpuFf^ Bf 
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(43) 
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(48) 



— _ Kii i) Fi'^ Bl j^KinFf^ B\ +_ KiiiFf^B\ 

These window functions were used and compared to the top hat window functions in [Rocha et al.| ( |20I0| l using the XFaster bandpower 
spectra and the Offset Lognormal Bandpower likelihood. However as the XFaster likelihood is estimated for each I, a comparison of 
observed to theoretical power spectrum does not make use of such windows. Instead the raw pseudo-Cf of the observations, the kernels, and 
transfer or filter functions are all that is required for such comparison (see section[43). 



3.3 The algorithm 

The power spectrum is estimated by the following procedure: 

• Generate Monte Carlo simulations of time-ordered data (TOD) for both signal and noise. The noise must have the same characteristics 
as the observed data, and in practice must be determined from the observed data. The simulated signal, on the other hand, can be almost 
anything, as it is simply a tracer of the effects of time and spatial domain filtering in the process, and used to calculate the transfer function 
Fi. In practice, it is convenient to use an approximate model of the CMB to generate the signal. 

• Make maps of the TOD using the same mapmaking code as used for the observations. 

• Estimate the pseudo-C, spectra and the spherical harmonic coefficients from the signal only maps to get the transfer function F^. 

• Estimate the pseudo-C, spectra from the noise only maps to compute the noise bias {Ni). The pseudo-spectra can be computed with 
anafast of HEALPix package ( [Gorski et al.| ( |2005] >) when the masks are the same for temperature and polarization, otherwise we use a 
specific code from the suite of XFaster modules. 
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• Iterate Equations [8] and [TO] to obtain an estimate of gi,. In tiie diagonal, isotropic approximation of XFaster, tlie computational cost of 
the iterative estimator is very small compared with that of the TOD generation and map-making stages. 

• The iterative estimator yields the Fisher information matrix automatically. An estimate of the band power covariance is given by , 
therefore we automatically get the uncertainty on the estimator. 

Large ensembles of signal+noise simulations are not required to estimate the band power covariance matrix as in the Master procedure, 
cutting the cost of Monte Carlo simulations by 1/3. Furthermore the covariance is not biased by an assumed model (which at very least 
requires the Master procedure to be run twice to be close to unbiased errors). 

As mentioned in Section [3TT] XFaster can estimate both auto-spectra and cross-spectra jointly, using the full covariance of the aim, 
via a multiple-map analysis. 

3.3.1 Computational scaling 

The overall scaling for XFaster without accounting for the signal and noise Monte Carlo simulations should go as Imaxirimaps x ripoi)^ 
for the internal Fisher calculation, where Upoi is either 1 or 3, with a further scaling of (jibi^s)^ for the outer iteration step. 

Currently the code is not optimised for speed. It could be sped up substantially by parallelizing the Fisher computation and would then 
scale linearly with number of processors. 

Approximate times for a single CPU for Phase2, with CTP binning and using 30 Fisher iterations are as follows: 

• aim and Ct from the 100 Monte Carlo simulated maps: ~ 8 hours (~ 5 minutes each) 

• kernel: 30 minutes 

• transfer function: 30 minutes (less if one relaxes the binning) 

• power spectrum: 1 hour 

• average mode to check for possible bias: 15-20 minutes 



4 RESULTS 

4.1 Results: Power spectrum 

We estimated transfer (filter) functions, kernels, and the power spectrum for the observed map described in Section[2] We also computed the 
power spectrum for the average of the signal-l-noise simulated maps. This average mode run checks for possible biases of the power spectrum 
estimator itself. In principle, the estimator if unbiased should follow closely the input signal Ct model used to generate the signal simulations. 

Figure [5] shows the power spectra estimated for the observed map for Phase 1 and Phase 2 and their la error bars. The power spectra 
recover accurately the input power spectrum in the middle range of multipoles, 30 ^ ^ ^ 1000. At high-^ ( ^ ^ 1000), they are impacted by 
the noise. At low-^, as the large scale structure of the observed map is a WMAP constrained realization, the estimated power spectrum is not 
necessarily consistent with, and exhibits a dispersion around, the best-fit spectrum. Comparing the power spectra for the diverse phases we 
conclude that the cut-sky anisotropic noise case (Phase lb) exhibits slightly greater uncertainties and slightly larger multipole to multipole 
variations at low-^ than the full-sky, isotropic noise case (Phase la). This is expected, as the cut-sky will induce correlations at low-^. The 
kernels correct these correlations; however, there is still a small residual dispersion. On the other hand the anisotropic noise will enhance 
the overall white noise level increasing the power spectrum uncertainty. For Phase 2 the dispersion of the power spectrum and uncertainties 
at low-^ are enhanced due to the residuals of correlated 1// noise. As in Figure |5] the power spectrum for Phase 1 is estimated from maps 
generated with a smaller number of detectors (four) than those for Phase 2 (twelve), and is therefore noisier by a factor of the order ^ . 
Therefore the error bars of the Phase 2 power spectrum are smaller than those of Phase 1 . The right thing to do though is to compare the 
power spectra estimated for the same number of detectors. In Figure[6]we compare the power spectra for Phase la (middle plot) and Phase 2 
(right hand side), both estimated on maps generated with all twelve detectors. As expected, the uncertainty for Phase la is now smaller than 
that of Phase 2. 

The power spectra for Phase 2 for both the symmetric and asymmetric beams are highly consistent. Therefore we conclude that the 
beam asymmetry is reasonably well-handled by XFaster (see Section [4.1.1| for more details). 

Figure [7] shows the transfer functions for Phase 1 and Phase 2. Except for Phase 2b, all are close to 1. This is expected, as we do not 
pre-filter the TOD and the only effect at low-^ is that due to the mapmaking step and to the limited number of Monte Carlos of the signal 
maps available. This means there is remaining sample variance on the large-scales although not very significant. We included the transfer 
function estimates because they are an integral part of the method and in real life they will not be equal to 1 . However to investigate and show 
that the significance of these small departures from 1 are not significant we estimated the power spectrum with transfer function=l as plotted 
in Figure [6] We also highlighted the fact that the transfer functions are estimated consistently between polarization types (i.e. taking into 
account cross-correlation between polarization modes in any realization). Note that in this work we did not take into account the remaining 
MC error in the transfer function in the final estimate of the power spectrum since any production run used in the real case will include 
many more realizations than used in this work (they could be easily included by adding to the final Fisher matrix if needed). However, for 
the asymmetric beam case the transfer function exhibits an upturn at high-^. This upturn tries to correct the mismatch between the 'real' 
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asymmetric beam and our assumed symmetric beam (as discussed in Section [4.1.1^ . As the input BB power spectrum model of the signal 
Monte Carlo simulations for Phase 2 is set to zero, we cannot constrain the BB transfer function. The transfer function obtained reflects the 
inadequacy of the input model and hence is close to zero. When estimating the power spectrum we replace the BB transfer function by the 
EE transfer function. 

Figure [S] shows the power spectrum estimated for the average mode run. Whereas Figure |9] shows the difference between the power 
spectrum obtained for the average mode run (with error bars) and the fiducial model, Ct, used as input, for phase2, asymmetric beam case. 
The stepwise decreases of the amplitude of the error bars are caused by the changes of the d's bins size. These plots show that the power 
spectrum follows the input signal Ct, confirming that XFaster is an unbiased estimator. These results were obtained using 100 Monte 
Carlo simulations. Considering 500 simulations for phase la the small departures of the transfer function from 1 reduces slightly. Due to 
computational constraints this was not feasible for phase2. However increases in computational capability over time make it now possible to 
generate thousands of simulations. Results for 143GHz channel will be presented in our upcoming papers [Rocha et ar] ( |2010[ l; |Ashdown et| 
[aTl pOTOl l. 

Figure [TO|shows a variation of the above plots in which the BB power spectrum replaces the noise Monte Carlo simulations. These plots 
show that whenever the noise Monte Carlo simulations exhibit issues in the sense that they do not reflect accurately the noise characteristics 
of the observations, replacing them by the BB power spectrum is an adequate procedure, as for Planck at 70 GHz the BB power spectrum is 
mostly dominated by noise. 

The power spectrum estimated for the observed map has been compared with those from several other methods jAshdown et al.|p010^ ). 
To further assess the power spectrum estimator we propagated this analysis to cosmological parameter estimation using the new XFaster 
likelihood as described in Section l431 



4.1.1 Symmetric and asymmetric beams 

To study the impact of beam asymmetries on power spectrum and cosmological parameter estimation, we took a minimal, non-informative 
approach. When computing the power spectrum for the observed maps convolved with symmetric and asymmetric beams we assumed a 
FWHM = 14' symmetric beam for both cases. 

We started by investigating the effect of this assumption on the shape of the transfer function, Fe, in Equations [6] and [7] or Ft in 
Equations[37][38] and[39] In effect we compute not only Fi, but also FpSbe, where Sbi is the correction to the beam transfer function Bp. This 
correction arises from assuming a symmetric beam when estimating the power spectrum of an observed map that in fact has been convolved 



with an asymmetric beam. Call this function the generalized transfer function, {BF)t. As mentioned in Section 4.1 since we do not pre-filter 
the TOD the only effect at low-^ is that due to the mapmaking step. Therefore the transfer function should be very close to unity, particularly 
for the symmetric beam. Figurejisjshows {BF)e for both cases. For the symmetric case this function is very close to 1 as expected, with tiny 
oscillations around 1 at low-£ consistent with our expectations. Apart from the transfer function for the BB mode for phase2 which is close 
to zero. Since the input BB power spectrum model of the signal Monte Carlo simulations for phase2 is set to zero we cannot constrain the 
transfer function for the BB mode. We resort to using the transfer function for the EE mode to estimate the BB power spectrum instead. 

The symmetric and asymmetric [BF)i differ. In the asymmetric case, the function exhibits an upturn at high i. This upturn tries to 
correct the mismatch between the 'real' asymmetric beam and our assumption. Hence the resulting power spectrum is pretty consistent 
for both cases, as shown in Figure [IT] Considering the input symmetric beam with FWHM = 14' and the estimated FWHM~ 13' for the 
asymmetric beam, this effect should be approximately (14/13)^, although anisotropic pixel filtering will add an extra component of aliasing 
in the maps at high-^. 

Although the power spectra look consistent, the parameter estimation shows departures of the order of (j/2 for some of the parameters, 
as shown in Section [43] To investigate this further we enhanced the previous plot into Figure [12] Though there is very good agreement 
between the two power spectra, there is still a slight bias for the asymmetric beam case. This bias is consistent with the small differences in 
the estimated parameters, in particular for n^, and log[10^''^s] (see Section ' 



4.3 I 



Next we study the impact of using two different sets of Monte Carlo simulations on the estimation of the power spectrum for the 
asymmetric beam case. Phase 2b. One set is the Monte Carlo simulations for the symmetric beam (Phase 2a) and the other the correct Monte 
Carlo simulations for the asymmetric beam (Phase 2b) study. On the right hand side of Figure[TT]we plot the power spectrum estimated using 
both sets of simulations. The power spectrum of the observed map of Phase 2b (convolved with an asymmetric beam), estimated using the 
Monte Carlo simulations for Phase 2a (convolved with the symmetric beam), is biased high at high-^, as expected. 



4.2 Results: Likelihood 

Following Section 1^2] we use one dimensional slices as an approximation to investigate the Non-Gaussianity of the likelihood, sampling in 
each qb direction around the maximum likelihood solution . This approximation is adequate if the band powers are not heavily correlated. 
Note that the likelihood slices are estimated along the q^ band power deviations and not along the qi power deviations for each multipole I, 
and hence are affected by the binning procedure. These slices are plotted in Figures [T4]|15[|16||17| 

Figures[T4]and[T5] compare the XFaster likelihood to four other likelihood approximations, Gaussian, Lognormal, Offset Lognormal, 
and Equal Variance ( [Bond Jaffe & KnoxH2000) ; |Rocha et al.| ( ]2010^ ). A thorough account of these likelihood approximations is given 
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Figure 5. Power spectrum estimated with XFaster and Icr error bars. Top row — Pliase la (left hand side) and Phase lb (right hand side) for map generated 
with a quadruplet of detectors. Bottom row — Phase 2a symmetric beam (left hand side) and Phase 2b asymmetric beam (right hand side) for map generated 
with all twelve detectors,. The plot displays the estimated power spectrum (blue) of the observed map, overplotted with the fiducial model used as input in 
Phase 2 signal simulations, first year WMAP+CBl+ACBAR best fit model (black). 




Figure 6. Power spectrum estimated with XFaster and Icr error bars, considering transfer function=l for Phase la, for map generated with a quadruplet of 
detectors (left hand side), with all twelve detectors (middle) and for Phase 2a, map generated with all twelve detectors and convolved with a symmetric beam 
(right hand side). This plot displays the estimated power spectrum (blue) of the observed map, overplotted with the Ci fiducial model used as input in Phase 2 
signal simulations, first year WMAP+CBI+ACBAR best fit model (black). 
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Figure 7. Transfer Function: Top row - for Phase la (left hand side) and Phase lb (right hand side); Bottom row - for Phase 2a symmetric beam (left hand 
side) and Phase 2b asymmetric beam (right hand side). 



elsewhere (see for instance Rocha et al. ( 2010 1), here we give a brief account of their definitions. In what follows C means the measured or 
observed quantity C. 



j2000f ), is a likelihood that is Gaussian in the Ce i.e. 



The Gaussian approximation, i Bond Jaffe & Knox j 
P(C'|C) oc exp |-^(C - C)^5"'(C' - C) I 

where C is a vector of Ci values (and similarly C) and S"^ is the inverse signal covariance matrix. 
The Offset Lognormal likelihood, ( [Bond Jaffe & KnoxlpOOO) ), is given by: 

Pln{C\C) oc exp - zfM{i - z)| , 

where zi = \n{Ce + xe) and the matrix M is related to the inverse covariance matrix by 

Mw = {Ce + Xi)Su}{C'e + x'l) 

(The offset factors xi are simply a function of the noise and beam of the experiment.) 
The Equal Variance likelihood, jBond Jaffe & Knox]pOOO| l) is given by: 



InL 
with 

2 = In (gt + 



{l-{z~z))\ 



and 



G 



[e (1 — (Tz)] ^ with az = 



(49) 



(50) 



(51) 



(52) 



(53) 



(54) 
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Figure 8. Power spectrum obtained with XFaster in average mode, meaning tliat we replaced the observed map by the average of the signal+noise simulated 
maps (blue). Top row — Phase la (left hand side) and Phase lb (right hand side) for map generated with a quadruplet of detectors. Bottom row — Phase 2a 
symmetric beam (left hand side) and Phase 2b asymmetric beam (right hand side) for map generated with all twelve detectors, overplotted with the Ci fiducial 
model used as input in our Phase 1 signal simulations, first year WMAP best fit model (black) for Phase 1 and first year WMAP+CBI+ACBAR best fit model 
(black) for Phase 2. It serves the purpose of checking for possible biases of the power spectrum estimator — in principle the power spectrum estimated in 
average mode should follow closely the input signal Ci model used to generate the signal simulations (black). 



The noise offset is estimated using the equation of the maximum likelihood solution for the qb replacing the observed map with the 
average of the noise Monte Carlo simulation power spectra {Nej. 

We use the values of the power spectrum and Fisher errors estimated with XFaster for Phase 2b (map generated with all twelve 
detectors and convolved with the asymmetric beam). As we are primarily interested in the shape of the likelihoods not on the actual value 
of the peaks, the comparison of the slices is done assuming they all peak at the same value, ie we use the band power specta estimated 
with XFaster and the functional shape of the other high-1 likelihoods. This is the same to say that we first compute the band power spectra 
and apply the functional forms whereas XFaster likelihood slices are computed when estimating the band-power spectra. Figure [74| shows 
temperature slices of the XFaster joint temperature and polarization likelihood. Figure [T5] shows slices for TT (first column), EE (second 
column), and TE (third column) (more precisely — In L). At the lowest multipoles the approximations differ, but as we move towards higher £ 
all but the Gaussian likelihood converge to the same functional form. For EE, however, the approximations differ noticeably for the Gaussian 
and Lognormal likelihood approximations (at I ~ 10, for instance). 

We further compare the XFaster likelihood to the exact likelihoods at low multipoles. The purpose is twofold. On one hand we want 
to validate the XFaster approximation, on the other hand we want to determine the £ range at which the approximations used for the high-£ 
power spectrum estimator breaks down. 

Figure [16] shows XFaster likelihood slices for the Phase 2 binned power spectra and the exact full-sky likelihood estimated for 
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Figure 9. Difference of the power spectrum obtained with XFaster in average mode, (meaning that we replaced the observed map by the average of the 
signal+noise simulated maps) and the Ci fiducial model used as input (first year WMAP+CBI+ACBAR best fit model), for phase!, asymmetric beam case. 
The stepwise decreases of the amplitude of the error bai's are caused by the changes of the bins size. 



Phase la. At low-^ the correlations induced by the cut-sky widen the XFaster likelihood, while the binning effect at high-£ results in a 
narrower distribution (when compared to the full-sky exact likelihood). Both effects are given by 

"=\/(2£-fi?/..,Al^^' ^^^^ 

where fsky is the fraction of the sky observed and is the width of the multipole bins. 

To make a direct comparison of the XFaster likelihood with the exact likelihood at low-^ estimated in the same cut-sky map (Phase 2), 
which accounts for the correlations induced by the sky cut, Figure[T7]shows the XFaster likelihood versus the pixel-based likelihood slices 
for temperature alone (Phase 2b, asymmetric beam case). The pixel based likelihood, BFLike, is a brute force likelihood evaluation of the 
multivariate Gaussian in pixel domain for a low-resolution map. The low-£ dataset of the CTP Phase 2 simulations was generated directly at 
A'^sidc ~ 16. In computing the slices, we conditioned on the remaining TT multipoles, Cj,"^ with £' / I, and on all multipoles of the TE 
and EE spectra (for details see |Rocha et al.| ( [2010| ). As for this case the BFLike estimates its own peak by computing a brute force pixel 
based likelihood on a downgraded map, we plot in Figure [TT] both likelihoods with its own peaks locations (left hand side) and assuming 
they peak at the same value (right hand side), ie after dividing both distributions by their peaks values. The plot on the left hand side might 
be misleading as the width of both distributions depend on their peaks locations. As we are mostly interested in comparing the shape of the 
likelihoods we will pay particular attention to the plot in the right hand side of Figure [TTl 

At ^ = 32 the agreement of both likelihoods is already quite remarkable, suggesting that a transition between high-^ and low-^ estimators 
around ^trans — 30^0 may be appropriate for this dataset. Hence a Planck hybrid likelihood built out of these two likelihoods (namely 
PiXFaster, see eg |Rocha et al.| ( [2010| l), with a transition range around i ~ 30^0 should be a viable hybridization scheme. 



4.3 Results: Cosmological parameters 

To compare the theoretical power spectrum with the observed power spectrum and estimate parameters, we need an operator to extract 
theoretical bandpowers from model power spectra Cj . In Rocha et al. |(|2010l we considered two types of windows, a top hat window per 



bin, and the appropriate Fisher- weighted window or XFaster band power window function. These window functions have been used in 
association with the Offset Lognormal Bandpower likelihood ( [Rocha et al.| ( |2010) ). The XFaster likelihood is estimated multipole by 
multipole ie for each I, hence no window function is required. In this mode XFaster can go straight from the map (via its raw pseudo-C^) 
to parameter estimation, bypassing the band power spectrum estimation step. 



We implemented the XFaster likelihood in a (suitably-modified) version of the publicly-available software CosmoMC| [ jLewis and 
[Bridle |f2002| l)) for cosmological parameter Markov Chain Monte Carlo estimations. The XFaster likelihood code computes the likelihood 
of a model passed to it by CosmoMC. There is no need for window functions or the band power spectrum itself. The inputs are the raw pseudo- 
Ct of the observations plus the kernel and transfer function required by XFaster to relate the cut-sky pseudo-C^ to the full-sky Ce. 



■* http://cosmologist.info/cosmomc/ 
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Figure 10. Power spectrum estimated with XFaster and Icr error bars for Pliase 2, for map generated witli all twelve detectors, for symmetric (top row) and 
asymmetric (bottom row) beams. Left hand side plot displays the estimated power spectrum (blue) if the observed map, overplotted is the Ci fiducial model 
used as input in our signal simulations, first year WMAP+CBI+ACBAR best fit model (black); Right hand side plot displays the estimated power spectrum 
with noise Monte Carlo simulations replaced by the BB power spectrum. 



Figures [18] through |20| show results for a simulation using Phase 2a (symmetric beams) and Phase 2b (asymmetric beams) data. The 
parameters considered are the baryon, cold dark matter and cosmological constant densities, u)^, = Q,^,h^ and cj^, = Q,ch^ and cja = ^a^^ 
respectively, the ratio of the sound horizon to the angular diameter distance at decoupling, 6s, the scalar spectral index n^, the overall 
normalization of the spectrum log[10^''y4] at fc = 0.05 Mpc~^ {As), the optical depth to reionization r, the age of the universe, the Hubble 
constant Ho, and the reionization redshift Zr^.- 

Figure[T9]shows parameters estimated for the average power spectrum of the signal+noise Monte Carlo simulations and for the observed 
power spectrum (i.e., the power spectrum estimated for the observed map). The parameters for the average simulated data recover the true 
input parameters, while those for the observed map shift from the input values, particularly for As- As mentioned before, the observed 
map is a WMAP-constrained realization, i.e., it uses the aim with phases measured by WMAP up to ^ = 70 to reproduce the large-scale 
structure observed by WMAP, and a best-fit model to the WMAP observations for I > 70. The WMAP best-fit parameters are obtained 
with considerable marginalization of the low-£ points by foregrounds. They are therefore unaffected by the low-^ anomalies. This means that 
unless we do such analysis too, we would not expect our observed realization to agree with the WMAP best fit model. This is clearly shown 
in Figure[T9] On the other hand, as Monte Carlo simulations are realizations of the WMAP best-fit model, one should expect no systematic 
bias from the ensemble of simulations, as confirmed in Figure [TS] 

As discussed in section 3, XFaster assumes that the noise is white (uncorrelated), i.e., that the noise covariance matrix is diagonal. 
Also, the XFaster likelihood is estimated multipole by multipole, hence to estimate the transfer function properly requires a larger number 
of Monte Carlo simulations to beat down the correlations between multipoles introduced by, e.g., sky cuts required for foreground removal. 
These simulations include both correlated noise and a sky cut. To assess whether the low-i? inadequacy of the likelihood is indeed the cause 
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Figure 11. Power spectrum estimated witli XFaster and la error bars for Piiase 2. Tiie leftiiand plot displays the power spectrum for a map generated 
with all twelve detectors, for both symmetric (red) and asymmetric (blue) beam cases. For both runs a symmetric beam with FWHM = 14' is assumed. The 
resulting power spectra are highly consistent, the compensation is achieved via the generalized transfer function (BF)^. The righthand plot displays the power 
spectrum estimated for map generated with a quadruplet of detectors, with two different sets of Monte Carlo simulations. In one case we use the Monte Carlo 
simulations convolved with the symmetric beam (blue); in the other we use the correct Monte Carlo simulations convolved with the asymmetric beam (red). 
Making use of the Monte Carlo simulations for the syimnetric case gives rise to a bias high at high-i?. 



of the parameter offsets seen in Figure[T9] we recalculated parameters, this time fixing r to the input model value in the simulations. Since r 
and only r is constrained almost entirely by the < 30 data, by fixing r we mimic the effect of using a likelihood evaluator that takes full 
cognizance of correlations in the noise and between multipoles at low I. The results are shown in Figure[T9] The estimated parameters shift 
towards the input parameter values, recovering those estimated for the average case in agreement with our postulated hypothesis. 

Table[T|compares parameter constraints for the symmetric beam case for the ensemble average power spectrum of the Monte Carlo sim- 
ulations run (avg) and the observed power spectrum run (pse) without r fixed. The parameter constraints tabulated are from the marginalised 
distributions. The parameter constraints for the observed power spectrum run (pse) when r is fixed to the fiducial input value are indistin- 
guishable from those derived from the average run, hence very close to the input parameter values. Therefore we do not include them in 
Tabled 

Figure [20lshows constraints from symmetric and asymmetric beam case for the average power spectrum of the Monte Carlo simulations 
and the actual observed estimated power spectrum without fixing r. Most of the parameters for both cases are consistent with each other. 
Investigating the plot for the average mode, we see deviations of the order of a/2 for Qch^, os, Us and Hq. There is an obvious degeneracy 
between erg and n^. For the observed case these deviations are noticeable mostly in As and cts. Once again these parameters are degenerate. 

The overall agreement in the parameter constraints from both symmetric and asymmetric beam cases is quite impressive. This reflects 
the adequacy of our procedure when dealing with beam asymmetries. 



5 CONCLUSIONS 

The XFaster power spectrum estimator is fully adequate to estimate the power spectrum of Planck data in the high-^ regime. It also 
performs well at moderately low multipoles, as long as the low-^ polarization and temperature power is properly accounted for, e.g., by adding 
an adequate low-^ likelihood ingredient. Our minimal non-informative approach enables us to recover most input parameters regardless of 
the asymmetry of the beam. 
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Figure 12. Power spectrum estimated witli XFaster and la error bars for Pliase 2, for map generated witii all twelve detectors, for both symmetric (red) and 
asymmetric (blue) beam cases. For both runs a symmetric beam with FWHM = 14' is assumed, the resulting power spectrum is consistent with each other, the 
mismatch is corrected via the generalized transfer function, {BF)i. There is still a slight bias for the asymmetric case in agreement with the small differences 
of the estimated cosmological parameters (see section [43) . 
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Figure 14. Likelihoods for Pliase 2b (12-detector map convolved with an asymmetric beam): Likelihood functions are sampled in each bandpower direction 
while fixing the other bands at the maximum likelihood values. The black (dotted) curve is the Gaussian approximation given by the Fisher matrix. The blue 
(dashed) curve is the offset lognormal approximation using the noise . The magenta (dash-dotted) curve is the equal variance approximation. The red 
(dashed) curve is lognormal distribution. The black (solid) curve is the XFaster likelihood estimated for temperature and polarization. The numbers in the 
right upper corner indicate the multipole I or ^effective of the binned multipoles. As for this set of multipoles h.i = 1 these numbers are the single multipole I 
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Figure 15. Likelihoods for Pliase 2b (12-detector map convolved with an asymmetric beam): Likelihood functions are sampled in each bandpower direction 
while fixing the other bands at the maximum likeHhood values. The black (dotted) curve is the Gaussian approximation given by the Fisher matrix. The blue 
(dashed) curve is the offset lognormal approximation using the noise . The magenta (dash-dotted) curve is the equal variance approximation. The red 
(dashed) curve is lognormal distribution. The black (solid) curve is the XFaster likelihood estimated for temperature and polarization, for TT (first column), 
EE (second column) and TE (third column); for = 5 (first row), £ = 7 (second row), £ = 10 for TT and TE and for bin with t in [10,13] for EE mode (third 
row). 




Figure 16. Likelihoods for Phase 2a (12-detector map convolved with a symmetric beam) for TT mode. XFaster likelihood (Phase 2 - cut sky) (blue solid 
line) vs Exact Full Sky likelihood (Phase la - full sky) (red solid line). At low-£ the XFaster likelihood is wider due to the correlations induced by the 
cut-sky while at high-£ XFaster likelihood is narrower due to the binning effect; The numbers in the right upper corner indicate the multipole I or the bin 
number. Up to 10 the bin number is the single multipole £ as A£ = 1, whereas bin = 50 corresponds to I in [61,62] and bin = 100 corresponds to £ in 
[257,263]. 
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Figure 17. Likelihoods for Pliase 2a (12-detector map convolved with a symmetric beam) for TT mode. Left hand side plot: XFaster likelihood slices (blue 
solid line) versus BFLike, Pixel based likelihood slices (red solid line), both curves are estimated on the cut-sky map of Phase 2. This plot can be misleading 
as the width of both distributions depend on their peaks locations. Right hand side plot: XFaster likelihood slices (blue solid line) versus BFLike, Pixel 
based likelihood slices (red solid line), both curves are estimated on the cut-sky map of Phase 2, assuming they both peak at the same value, ie after dividing 
both distributions by their peaks values. The agreement is already apparent at I as low as £ = 16, 32. The numbers in the right upper comer indicate the 
multipole I or ^effective of the binned multipoles. As for this set of multipoles A£ = 1 these numbers are the single multipole I 



Table 1. Parameter constraints from Planck 70 GHz data. Phase 2a, for map generated with all twelve detectors and convolved with a symmetric beam, 
using XFaster power spectrum and likelihood estimator for the average of the signal-l-noise Monte Carlo simulations power spectrum (avg) and the actual 
observed power spectrum (pse) without r fixed. As mentioned in the text, the parameter constraints for the actual power spectrum run (pse) when r is fixed 
to the fiducial input value are indistinguishable from those derived from the average run, hence very close to the input parameter values. Parameter constraints 
displayed here were obtained from the marginalised distributions. 
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Figure 18. Parameter constraints from Phase 2a (12-detector map convolved with a symmetric beam). The 1 -dimensional marginalized posteriors are from 
XFaster likelihood (with inclusion of all modes, TT, EE, BB, TE): for ensemble average of signal+noise Monte Carlo simulations, i.e., average run (solid 
black lines), for several single signal+noise Monte Carlo simulations, and for values of the fiducial best fit input parameters (black vertical lines). The 
parameters for the average power spectrum recover the true input parameters. Furthermore this plot shows that there is no systematic bias for each Monte 
Carlo simulation. This is to be expected as the Monte Carlo simulations are realizations of the WMAP best fit model. 
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Figure 19. Parameter constraints from Phase 2a (12-detector map convolved wfith a symmetric beam). The 1 -dimensional marginalized posteriors are from 
XFaster likelihood (with inclusion of all modes, TT, EE, BB, TE): for ensemble average of signal+noise Monte Carlo simulations, i.e., average run (solid 
black lines), for single observed data (dashed blue lines) and for the run with t = Tfiducial using XFaster likelihood for ^ > 30 (solid red lines), 
with overplotted values of the fiducial best fit input parameters (black vertical lines). The parameters for the average power spectrum recover the true input 
parameters. However the parameters for the observed map shift from the input parameters, particularly so for the parameter Amplitude As - As the large scale 
structure of the observed map is a WMAP constrained realization we do not expect the estimated parameters to agree with WMAP best fit parameters. Fixing t 
to the input value regularizes the amplitude in the likelihood runs, now the estimated parameters of the observed map shift towards the input parameter values 
recovering those estimated for the ensemble average of the Monte Carlo simulations. 
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Figure 20. Parameter constraints from Phase 2a (12-detector map). The 1 -dimensional marginahzed posteriors are from XFaster Hkelihood (with inclusion 
of all modes, TT, EE, BB, TE): for the averaged spectra (left hand side) and for single realization spectra (observed map) without fixing t (right hand side), 
asymmetric beam (solid blue lines) and symmetric beam (dashed red lines) cases. 
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